home *** CD-ROM | disk | FTP | other *** search
/ Almathera Ten Pack 2: CDPD 1 / Almathera Ten on Ten - Disc 2: CDPD 1.iso / pd / 201-225 / 214 / mandelvroom / src / mandffp.c < prev    next >
C/C++ Source or Header  |  1995-03-13  |  9KB  |  426 lines
  1. /*
  2.  * MandelVroom 2.0
  3.  *
  4.  * (c) Copyright 1987,1989  Kevin L. Clague, San Jose, CA
  5.  *
  6.  * All rights reserved.
  7.  *
  8.  * Permission is hereby granted to distribute this program's source
  9.  * executable, and documentation for non-comercial purposes, so long as the
  10.  * copyright notices are not removed from the sources, executable or
  11.  * documentation.  This program may not be distributed for a profit without
  12.  * the express written consent of the author Kevin L. Clague.
  13.  *
  14.  * This program is not in the public domain.
  15.  *
  16.  * Fred Fish is expressly granted permission to distribute this program's
  17.  * source and executable as part of the "Fred Fish freely redistributable
  18.  * Amiga software library."
  19.  *
  20.  * Permission is expressly granted for this program and it's source to be
  21.  * distributed as part of the Amicus Amiga software disks, and the
  22.  * First Amiga User Group's Hot Mix disks.
  23.  *
  24.  * contents:  this file contains the routines to calculate Mandelbrot and
  25.  * Julia pictures in Amiga Fast Floating Point.
  26.  */
  27.  
  28. #include "mandp.h"
  29.  
  30. extern SHORT MaxOrbit;
  31.  
  32. static LONG MaxIter;
  33.  
  34. /*
  35.  * Fast Floating Point Mandelbrot Generator
  36.  */
  37. MandelbrotFFP( Pict, StartX, StartY, GapX, GapY )
  38.   /*
  39.    * These parameters are IEEE 32 bit floating point numbers
  40.    */
  41.  
  42.   struct Picture *Pict;
  43.  
  44.   LONG  StartX, StartY, GapX, GapY;
  45. {
  46.   register int i, j, k;
  47.   register SHORT *CountPtr;
  48.  
  49.   register float curx, cury;
  50.  
  51.   float startx, starty, gapx, gapy;
  52.   float SPFieee();
  53.  
  54.   UBYTE MandFlag;
  55.  
  56.   /* Here we convert the IEEE parameters into FFP format */
  57.  
  58.   startx = SPFieee( StartX );
  59.   starty = SPFieee( StartY );
  60.   gapx   = SPFieee( GapX );
  61.   gapy   = SPFieee( GapY );
  62.  
  63.   MaxIter = Pict->MaxIteration;
  64.  
  65.   if (Pict->Flags & NO_RAM_GENERATE)
  66.     CountPtr = Pict->Counts;
  67.   else
  68.     CountPtr = Pict->Counts + Pict->CurLine*Pict->CountX;
  69.  
  70.   /* start in the upper left hand corner */
  71.  
  72.   cury = starty;
  73.  
  74.   /*
  75.    * for each row
  76.    */
  77.   for (i = Pict->CurLine; i < Pict->CountY; i++) {
  78.  
  79.     curx = startx;
  80.  
  81.     MandFlag = 0;
  82.  
  83.     if ( Pict->Flags & NO_RAM_GENERATE )
  84.       CountPtr = Pict->Counts;
  85.  
  86.     /*
  87.      * for each collumn
  88.      */
  89.  
  90.     for (j = 0; j < Pict->CountX; j++) {
  91.  
  92.      /* if we've hit a Mandelbrot point, then use the convergence detector,
  93.       *  to reduce compute time
  94.       */
  95.  
  96.       if (*CountPtr == 0) {
  97.  
  98.         if ( MandFlag ) {
  99.  
  100.           k = FFP_Trace_Height( 0.0, 0.0, curx, cury );
  101.  
  102.         } else {
  103.  
  104.           k = FFP_Height( 0.0, 0.0, curx, cury ); /* Normal calculator */
  105.         }
  106.  
  107.         MandFlag = k == Pict->MaxIteration;
  108.  
  109.         *CountPtr = k;        /* save it in ram for recoloring */
  110.       }
  111.       CountPtr++;
  112.  
  113.       curx += gapx;
  114.  
  115.       ChildPause( Pict );
  116.     }
  117.     cury += gapy;
  118.  
  119.     CheckEOL( Pict );
  120.   }
  121. } /* Mandelbrot */
  122.  
  123. /***************************************************************************
  124.  *
  125.  *           Julia Curve Generator that uses Amiga Fast Floating Point
  126.  *
  127.  **************************************************************************/
  128.  
  129. /*
  130.  * Fast Floating Point Julia Generator
  131.  */
  132. JuliaFFP( Pict, JuliaX, JuliaY, StartX, StartY, GapX, GapY )
  133.  
  134.   /*
  135.    * These parameters are IEEE 32 bit floating point numbers
  136.    */
  137.  
  138.   register struct Picture *Pict;
  139.   LONG  JuliaX, JuliaY;
  140.   LONG  StartX, StartY;
  141.   LONG  GapX,   GapY;
  142. {
  143.   register int i, j, k;
  144.   register float curx, cury;
  145.  
  146.   register SHORT *CountPtr;
  147.  
  148.   float juliax, juliay, gapx, gapy,startx;
  149.   float SPFieee();
  150.  
  151.   UBYTE ConvFlag;
  152.  
  153.   MaxIter = Pict->MaxIteration;
  154.  
  155.   if (Pict->Flags & NO_RAM_GENERATE)
  156.     CountPtr = Pict->Counts;
  157.   else
  158.     CountPtr = Pict->Counts + Pict->CurLine*Pict->CountX;
  159.  
  160.   /* Here we convert the IEEE parameters into FFP format */
  161.  
  162.   juliax = SPFieee( JuliaX );
  163.   juliay = SPFieee( JuliaY );
  164.  
  165.   gapx = SPFieee( GapX );
  166.   gapy = SPFieee( GapY );
  167.  
  168.   /* start in the upper left hand corner */
  169.  
  170.   startx = SPFieee(StartX);
  171.   cury   = SPFieee(StartY);
  172.  
  173.   /* move down to next not generated line */
  174.  
  175.   cury += Pict->CurLine*gapy;
  176.  
  177.   /*
  178.    * for each row
  179.    */
  180.   for (i = Pict->CurLine; i < Pict->CountY; i++) {
  181.  
  182.     curx = startx;
  183.     ConvFlag = 0;
  184.  
  185.     if ( Pict->Flags & NO_RAM_GENERATE )
  186.       CountPtr = Pict->Counts;
  187.  
  188.     /*
  189.      * for each collumn
  190.      */
  191.  
  192.     for (j = 0; j < Pict->CountX; j++) {
  193.  
  194.      /* if we've hit a Mandelbrot point, then use the convergence detector,
  195.       *  to reduce compute time
  196.       */
  197.  
  198.       if (*CountPtr == 0) {
  199.  
  200.         if (ConvFlag) {
  201.                                            /* try the convergence method */
  202.           k = FFP_Trace_Height( curx, cury, juliax, juliay);
  203.         } else {
  204.           k = FFP_Height( curx, cury, juliax, juliay);      /* Normal calc */
  205.         }
  206.  
  207.         ConvFlag = k == Pict->MaxIteration;
  208.  
  209.         *CountPtr = k;        /* save it in ram for recoloring */
  210.       }
  211.       CountPtr++;
  212.  
  213.       curx += gapx;
  214.  
  215.       ChildPause( Pict );
  216.     }
  217.     cury += gapy;
  218.  
  219.     CheckEOL( Pict );
  220.   }
  221. } /* JuliaFFP */
  222.  
  223. /*
  224.  * This function calculated the 'potential' of a given location
  225.  * in the complex plane.
  226.  */
  227. int
  228. FFP_Height( posx, posy, juliax, juliay )
  229.   float posx;    /* Location in screen coordinate space         */
  230.   float posy;    /* Location in screen coordinate space         */
  231.   float juliax;  /* Real part of location on complex plane      */
  232.   float juliay;  /* Imaginary part of location on complex plane */
  233. {
  234.   register float cura, curb;
  235.   register float cura2, curb2;
  236.   register LONG k;
  237.  
  238. #ifdef CHECK_TASK_STACK
  239.  
  240.   CheckStack();
  241.  
  242. #endif
  243.  
  244.   cura = cura2 = posx;
  245.   curb = curb2 = posy;
  246.  
  247.   cura2 *= cura2;
  248.   curb2 *= curb2;
  249.  
  250.   for (k = 0; k < MaxIter; k ++ ) {
  251.  
  252.     curb *= cura;                     /* b = 2 * a * b  + juliay    */
  253.     curb += curb + juliay;
  254.  
  255.     cura = cura2 - curb2 + juliax;    /* a = a2 - b2 + juliax       */
  256.  
  257.     cura2 = cura * cura;              /* a2 = a * a                 */
  258.     curb2 = curb * curb;              /* b2 = b * b                 */
  259.  
  260.     if (cura2+curb2 >= 16.0)          /* quit if a2+b2 is too big   */
  261.       return( k );
  262.   }
  263.   return( k );
  264. }
  265.  
  266. /*
  267.  * this is also a Mandelbrot potential calculator, that is tailored to
  268.  * identify points that are 'in' the Mandelbrot set.  Points outside the
  269.  * Mandelbrot set diverge.  Points inside converge.
  270.  * This code uses a trace table mechanism to detect convergence.
  271.  */
  272. int
  273. FFP_Trace_Height( posx, posy, curx, cury )
  274.  
  275.   float posx;    /* Real part of location on complex plane      */
  276.   float posy;    /* Imaginary part of location on complex plane */
  277.  
  278.   float curx;    /* Real part of location on complex plane      */
  279.   float cury;    /* Imaginary part of location on complex plane */
  280. {
  281.   LONG  k, l;
  282.  
  283.   LONG TraceSize = 32;
  284.  
  285.   float olda[32], oldb[32]; /* Convergence trace table */
  286.  
  287.   register float cura, curb, cura2, curb2;
  288.   register float *RealTrace, *ImagTrace;
  289.  
  290. #ifdef CHECK_TASK_STACK
  291.  
  292.   CheckStack();
  293.  
  294. #endif
  295.  
  296.   cura = cura2 = posx;
  297.   curb = curb2 = posy;
  298.  
  299.   cura2 *= cura2;
  300.   curb2 *= curb2;
  301.  
  302.   for (k = 0; k < MaxIter; k += TraceSize) {
  303.  
  304.     RealTrace = olda;
  305.     ImagTrace = oldb;
  306.  
  307.     for (l = 0; l < TraceSize; l++) {
  308.  
  309.       *(RealTrace++) = cura;
  310.       *(ImagTrace++) = curb;
  311.  
  312.       curb *= cura;
  313.       curb += curb + cury;
  314.  
  315.       cura = cura2 - curb2 + curx;
  316.  
  317.       cura2 = cura * cura;
  318.       curb2 = curb * curb;
  319.  
  320.       if (cura2+curb2 >= 16.0)
  321.         return( k + l );
  322.     }
  323.  
  324.     /* Scope out trace table for convergence */
  325.  
  326.     RealTrace = olda;
  327.     ImagTrace = oldb;
  328.  
  329.     for (l = 0; l < TraceSize; l++)
  330.  
  331.       if (cura == *(RealTrace++) && curb == *(ImagTrace++)) {
  332.         k += MaxIter;
  333.         return( MaxIter );
  334.       }
  335.   }
  336.   return( MaxIter );
  337. }
  338.  
  339. DrawOrbitFFP( Pict, creal, cimag, zreal, zimag )
  340.   register struct Picture *Pict;
  341.   LONG zreal, zimag, creal, cimag;
  342. {
  343.   register struct RastPort *Rp;
  344.   register float cura,  curb;
  345.            float cura2, curb2;
  346.            float Creal, Cimag;
  347.  
  348.   float x_scale, y_scale;
  349.   int x_center, y_center;
  350.   int width, height;
  351.   int x, y;
  352.   register int k;
  353.  
  354.   struct Window *Window;
  355.  
  356.   Window = OrbitWind;
  357.  
  358.   Rp = Window->RPort;
  359.  
  360.   width  = (Window->Width-Pict->LeftMarg-Pict->RightMarg);
  361.   height = (Window->Height-Pict->TopMarg-Pict->BotMarg);
  362.  
  363.   x_center = width/2 + Pict->LeftMarg;
  364.   y_center = height/2 + Pict->TopMarg;
  365.  
  366.   y_scale = x_scale = (float) height / 2.0;
  367. /*x_scale *= AspectRatio( Pict );*/
  368.  
  369.   if ( Pict->pNode.ln_Type == MANDPICT ) {
  370.  
  371.     Creal = SPFieee(creal);
  372.     Cimag = SPFieee(cimag);
  373.  
  374.     cura = cura2 = SPFieee(zreal);
  375.     curb = curb2 = SPFieee(zimag);
  376.  
  377.   } else {
  378.  
  379.     Creal = SPFieee(zreal);
  380.     Cimag = SPFieee(zimag);
  381.  
  382.     cura = cura2 = SPFieee(creal);
  383.     curb = curb2 = SPFieee(cimag);
  384.  
  385.   }
  386.  
  387.   cura2 *= cura2;
  388.   curb2 *= curb2;
  389.  
  390.   SetAPen( Rp, 0 );
  391.   RectFill( Rp, Pict->LeftMarg, Pict->TopMarg, width+Pict->LeftMarg-1, height+Pict->TopMarg-1);
  392.  
  393.   SetAPen( Rp, HIGHLIGHTPEN);
  394.  
  395.   for (k = 0; k < MaxOrbit; k++ ) {
  396.  
  397.     curb *= cura;
  398.     curb += curb + Cimag;
  399.  
  400.     cura  = cura2 - curb2 + Creal;
  401.  
  402.     cura2 = cura * cura;
  403.     curb2 = curb * curb;
  404.  
  405.     if (cura2+curb2 >= 16.0)
  406.       return( k );
  407.  
  408.     /* map real and imaginary parts into window coordinates */
  409.  
  410.     x = x_center + (int)(x_scale*cura);
  411.     y = y_center + (int)(y_scale*curb);
  412.  
  413.     if ( x >= Pict->LeftMarg && x < Window->Width-Pict->RightMarg &&
  414.          y >= Pict->TopMarg  && y < Window->Height-Pict->BotMarg ) {
  415.  
  416.       /* plot pixel location */
  417.       WritePixel( Rp, x, y);
  418.     }
  419.  
  420.   }
  421.   return( k );
  422. }
  423.  
  424.  
  425.  
  426.